1 

Final State of Gregory- Laflamme Instability 

Luis Lehner and Frans Pretorius 



o 

£3 ■ Abstract 

We describe the behavior of a perturbed 5-dimensional black string subject 

■ to the Gregory-Lafiamme instability. We show that the horizon evolves in a 
self-similar manner, where at any moment in the late-time development of 
the instability the horizon can be described as a sequence of 3-dimensional 

Q-i spherical black holes of varying size, joined by black string segments of simi- 

lar radius. As with the initial black string, each local string segment is itself 
unstable, and this fuels the self-similar cascade to (classically) arbitrarily 
small scales; in the process the horizon develops a fractal structure. In finite 
asymptotic time, the remaining string segments shrink to zero-size, yielding 
a naked singularity. Since no fine-tuning is required to excite the instability, 
this constitutes a generic violation of cosmic censorship. We further dis- 
cuss how this behavior is related to satellite formation in low-viscosity fluid 
streams subject to the Rayleigh-Plateau instability, and estimate the fractal 

■ dimension of the horizon prior to formation of the naked singularity. 
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1.1 Overview 



The ultimate fate of black holes subject to the Gregory-Lafiamme instability 
has been an open question for almost two decades. In this chapter we discuss 
the behavior of an unstable 5D black string and elucidate its final state. Our 
studies reveal that the instability unfolds in a self-similar fashion, where the 
horizon at any given time can be seen as thin strings connected by hyper- 
spherical black holes of different radii. As the evolution proceeds, pieces of 
the string shrink while others give rise to further spherical black holes, and 
consequently the horizon develops a fractal structure. At this stage its overall 
topology is still R x S 2 , with the fractal geometry arising along R and hav- 
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ing an estimated Hausdorff dimension d ~ 1.05. However, the ever-thinning 
string-regions eventually shrink to zero size, revealing a (massless) naked 
singularity. Consequently, this spacetime provides a generic counterexample 
to the cosmic censorship conjecture, albeit in 5 dimensions. While we re- 
strict to the 5D case for computational cost reasons, our observations are 
intuitively applicable to higher ones. 

To correctly capture the late-time nonlinear dynamics of the system re- 
quires numerical solution of the full Einstein equations. In this chapter, 
following a brief historical account (Sec. I1.2p . we describe details of our nu- 
merical implementation (Sec. 1 1 . 3|) as well as the behavior of the obtained 
solution (Sec. II. 4|) . We discuss some additional properties of the solution, 
including speculation on when quantum corrections are expected to become 
important, and future directions, in Sec. 11.51 



1.2 Background 

Gregory and Laflamme's observation ([Gregory and Laflammel . 1993) that 
linearized perturbations of D-dimensional black strings (D > 5) are unsta- 
ble for long wavele ngths is in stark contrast with the known behavior of 
black holes in 4D (jPricd . Il972l ). Based on the nature of the growing per- 
turbations of the black strings, together with entropy arguments, Gregory 
and Laflamme conjectured black strings would bifurcate, thus inducing a 
topology change in the horizon to yield localized S D ~ 2 black holes. How- 
ever, black hole bifurcation necessarily implies, at the classical level, a for- 
mation of a naked singu larity where the pinch-off occurs (see for example 
Hawking and Ellisl . Il973h . Assuming such behavior would be resolved by 



quantum gravity, the conjecture was taken as likely true for about a decade. 
In early 2000s, tension arose when Horowitz and Maeda proved that any 
bifurcation could only take place at infinite affine t ime along the genera- 
tors of the horizon that cross the bifurcation point ([Horowitz and Maedal . 



200 ll ). They dismissed this possibility as unlikely and conjectured the exis- 



tence of stationary, non-uniform black string solutions which would be the 
end point of the instability. Follow up works pres ented approxima te station- 
ary solutions found perturbatively ( Gubser . 20021 ) or numerically ( Wiseman . 



2003 ; Kudoh and Wiseman . 20041 ) . though these had less entropy than the 



uniform string and so could not be the end-point of the system. To these 
developments, interesting observat ions we r e ma de as to a possible rever- 
sal of this behavior for D > 1 3 (ISorkinl. liooih but that this dimension 
varies in boosted black strings ( Hovdebo and Mvers . 20061 ). Furthermore, 
it was pointed out that electrically charged black strings are more unstable 
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than magnetically charged ones (jSarbach and Lehnerl . 120051 ) , and arguments 
were presented for a conical s t ructure in the black string-black hole transi- 
tion (|Kq1 and Wisemanl . 120031 : Eol I2OO6I ). 

This flurry of activity not only hinted at the possibility of rich phe- 
nomenology awaiting in the dynamics of the system, but the need for a full 
analy sis to uncover it. A first attempt to do so was presented in lChoptuik et al 
20031 . This study revealed that the development of a black string, perturbed 
by a long wavelength periodic mode, progressed to a structure that could be 
described as a sequence of S 3 black holes joined by strings, though could not 
uncover the final fate as the code was unable to evolve the solution further. 
Detailed follow-up analysis of the results from this simulation showed that 
the affine time al ong the generators gre w faster than a simple exponential of 
asymptotic time ( Garfinkle et al. . 20051 ) , suggesting consistency of a possible 
pinch-off with the theorem presented in lHorowitz and Maedall200ll (see also 
Marolj [20051 1. 

Following these works, additional hints to the possib le end state came by 
analogy with fluid systems. The membrane paradigm ( Thorne et all Il986l ) 



first suggested that event horizon dynamics could be described, to lead- 
ing order, by the Navier-Stokes equations for a viscous fluid (albeit with 
some "unusual" properties, such as a negative bulk viscosity). More re- 
cently, descriptio ns of black brane dynamics in asymptotically 5D Anti de- 
Sitter spacetime ( Bhattacharvva et al. . 20081 ) . and an effective worldyolume 
theor y of black holes (the "black folds" paradigm, lEmparan et all 2009, 



2010l ). similarly describe the dynamics by Navier-Stokes equations. Cardoso 
and Dias noted the qualitative similarity between the dispersion relation of 
unstable modes of black strings to those of thin fluid streams subject to 
the Raleigh-Plateau instability 0, and pointed out that if the similarities 
persisted beyond linear development the black s tring solution could shar e 
the same fate as the fluid, which does pinch-off ( Cardoso and Dias . 20061 ). 
Moreover, in the latter system evolution to pinch-off can be accompanied by 
a phenomena known as satellite formation, where one or more generations 
of ever smaller spherical "beads" form in the thinning stream; the lower 
the viscosity of t he fluid, the more generations are observed (see for exam- 
ple lEggersl . Il997l ). In fact, as we d escribe here, this is qualitat ively what 
is seen to occur in the black string ( Lehner and Pretorius . 2010l ). The self- 



1 Though note that the relationship is not exact; in particular in the long wavelength A limit 
the growt h rate il in Raleigh-Plateau differs from Greg ory-Laflamme as f2 oc A - 1 / a vs 
fl oc A -1 (Cardoso and Dias, 2006; Camps ct al., 2010). However, as the results here show, 
the behavior of the system is determined by "intermediate" wavelengths, thus the infinitively 
long wavelength behavior is likely irrelevant as far as the late time behavior of the string is 
concerned. 
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similar nature of this solution makes it impossible to numerically follow the 
evolution to arbitrary small scales, though following the trends, together 
with the fact that there is no intrinsic length scale in the field equations 
of general relativity, suggests that infinitely many generations of satellites 
will form (at the classical level). Interestingly, the horizon- fluid analogues 
describe the horizon as a perfect fluid with a shear viscosity to entropy 
densi ty ratio rj/s of 1/4-7T, implying a lowe r visco sity than any "real-world" 
fluid dKovtun et all 1200.4 buchel and Lhl I2004T ). 



1.3 Numerical approach 

As mentioned, we relied on numerical simulations to reveal the solution of 
an unstable black string. In this section we describe particularly relevant 
aspects of the numerical implementation adopted. We note however that 
achieving a reliable implementation of Einstein's equations for any given 
class of problem require dealing with many subtle issues. It is beyond the 
scope of this chapter to go into the full details (for textbook introductions to 
numerical relativity see | Bona and Palenzuela-Luqud . 120051 : lAlcubierrd . 120081 : 
Baumgarte and Shapira 2O10l ): however, since there is as of yet no universal 
approach to all problems requiring numerical relativity, we find it useful to 
describe the method that has been successful in evolving the 5D black string 
spacetime, focusing on details particular to this problerrJl. 

In the class of "not-so-subtle" issues for successful solution of the ini- 
tial boundary value problem are: solving the constraint equations for con- 
sistent initial data; using a mathematically well-posed formulation of the 
field equations that furthermore does not admit exponentially growing so- 
lutions from truncation-error-seeded constraint violations; dealing with the 
geometric singularity inside the black string; and using a stable numeri- 
cal integration scheme that can adequately resolve all the relevant length 
scales. Th e initial data we use here is the same as in the earlier numeri- 
cal study (|Choptuik et all hooi ). and we describe the ansatz and solution 
method in Se c.ll.3.Tl For evolu tion, the present study uses the harmonic de 
osition (Pretoriusl . l2005bfl with constraint damping (| Gundlach et al 



composition (i rretorius, 
200.4 IPretoriusl . 1200581 1 



and is described in Sec. 11.3.21 We use the excision 
approach to remove the black string singularity from the computational do- 
main. This relies on having a detailed description of the apparent horizon 

2 For other examples of nu merical relativity exploring dynamical scenarios in higher 

dimensional settings see (Choptuik ct al., 2003; Sarbach and Lclmer, 2004; Guzman et al., 



(let 

120071; IWitek et all.i2010l. 120111; IShibata and Yoshind 1201(3 ; IQkawa et all 1201 ID 
3 For additional examples of its use in numerical relativity see, e.g. ( Garfinklc, 120021 ; 
ISzilagvi and Winicouil 120031 : iLindblom et all l2006t IPalenzuela et all I2009F 
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of the spacetime, which is also a key structure employed to analyze the dy- 
namics of the unstable string. In Sec. IL3.3l we describe basic properties of an 
apparent horizon and how it is found in the simulation. For our numerical 
scheme, we use finite difference techniques with adaptive mesh refinement 
(AMR). These methods will not be discussed here, but where appropriated 
references to more information will be given. 

"Subtle" issues of import to the 5D black string evolution will be discussed 
in the following sections as well, including the coordinates employed and an 
appropriate choice of constraint damping parameters. 



1.3.1 Initial Data 

We are interested in studying the black string's dynamics with respect to per- 
turbations in vacuum. Here, rather than adopting generic perturbations we 
concentrate on those that only break the symmetry along the extra "string" 
direction w (henceforth we will re fer to it as "string- like" direct i on). A s has 
been shown in the original work (j Gregory and Lanammel . Il99.4 Il994f ) , per- 



turbations in the S 2 (6, 4>) sector lead, at the linear level, to sub-dominant 
modes with respect to those in w. Thus restricting to such a case should not 
affect the generality of the observed behavior and, as an important added 
bonus, allows us to restrict to a numerical implementation with only 2 rele- 
vant spatial coordinates. 

In what follo ws we summariz e the a pproach adopted to define such data (for 
full details see Choptuik et al. . 20031 ). To obtain consistent initial data, we 



adopt a Cauchy (4 + 1) decomposition of Einstein equations. At an initial 
(t = 0) hypersurface, its intrinsic metric ("fij) and extrinsic curvature (Kij) 
provide suitable data if they satisfy the Hamiltonian and momentum con- 
straints: 

H = { ^R + K 2 - K i:j K ij = 0, (1.1) 
M i = D j (Kf-iK) = 0, (1.2) 

with (^R the Ricci scalar associated to jij, K = -y^Kij and Di the covariant 
derivative compatible with 7^. We only consider perturbations depending 
on (r,w), respecting the SO (3) symmetry; then the general intrinsic metric 
can be expressed as, 

^ds 2 = ~{ rr dr 2 + 2y TW drdw + jwwdw 2 + r 2 yndQ, 2 , (1-3) 

where dQ 2 is the unit 2-sphere metric. With these assumptions Mg = = 
so the constraints provide three equations for eight variables. Thus three un- 
knowns can be solved for provided the remainder information is prescribed. 
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To do so, notice that in ingoing Eddington-Finkelstein coordinates, the un- 
perturbed black string has 

2M . . 

7rr = H , 7ru , = , j ww = 1 , 7n = 1 , (1.4) 

r 

K„ = -2M (r + 3 M) a , K rw = , K ww = 1 , = 2M« ; (1.5) 

with a = y/r/(r + 2M) and M the mass per unit length. We thus adopt 
Irw = Iww — 1 = ^riu = ^wiu = and introduce the perturbation by defining 

7n = 1 + A sin (w^j e-^-rofM . (L6) 

Here, A is a parameter controlling the overall strength of the perturba- 
tion, q is an integer to defining its spatial frequency in the w direction, and 
ro and S r are parameters controlling the extent of the perturbation in the 
radial direction. For the results presented here A = 0.1, q = l,r D = 2.5 
and 5 r = 0.5. The remaining variables are obtained numerically by solv- 
ing the Hamiltonian and momentum constraints. To this end, we adopt a 
finite difference approximation of the constraints on a uniform grid with 
(r, w) £ [Rmin, Rmax] x [0> L] and deal with the constraints as follows. The 
Hamiltonian constraint provides an equation for 7rr which can be schemat- 
ically written as 

F\ d r "f rr + F2 "frrd-ww"frr ~i~ F3 "frrdw'Jrr ~i~ F4 (^luTrr) F§ (Trr) -^isTVr — 0) 

(1-7) 

with Fi (I = 1..6) being functions that do not depend on 7rr . This equa- 
tion is integrated outwards from an inner boundary chosen well inside the 
unperturbed black string horizon {R m in = M), and boundary data there is 
provided by the unperturbed value of 7rr . A second order accurate radial 
integration is adopted where the presence of the w— derivatives make the 
numerical problem a non-linear, cyclic (due to the ^-periodicity), tridiago- 
nal system for the unknowns jrrlnw ^^is s y s * em is solved using Newton's 
method and a cyclic tridiagonal linear algorithm. 

The momentum constraint along the r direction provides a first order 
equation for Kgg of the form, 

ddrKge + G 2 K ee + G 3 = , (1.8) 

with Gi (I = 1..3) being functions that do not depend on Kgg. This is a simple 
ODE equation which is integrated along each w = const, line, employing 
second-order finite differences, and boundary conditions at R max provided 
by the unperturbed black string solution. The momentum constraint along 
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the w direction provides an equation for K rr of the form 

H\d w K rr + H 2 K„ + #3 = 0, (1.9) 

with Hi (I = 1..3) being functions that do not depend on K rr . This is 
again a simple ODE that is integrated to second order accuracy along r = 
const, lines. Boundary conditions are specified along w = w m i n with the 
unperturbed values. 

The procedure above provides consistent data for our problem. We stress 
that it is by no means a general procedure to study consistent initial data 
for a black string problem, rather it is a straightforward way of setting 
up a particular perturbation of the black string to study its subsequent 
evolutior0. 

The method described above provides data for {g^ , Kij} (i, j = r, 9, (f>, w). 
However, data required for the harmonic formulation (discussed in the fol- 
lowing sub-section) consists of {g^ u , 9fiu,t}- Note also that for evolution 
we impose the harmonic condition with respect to Cartesian coordinates 
(x,y,z,w), thus we first transform (r,9,<fi) to (x,y,z) via the standard re- 
lations between 3D Cartesian and spherical polar coordinates. Data for 
can be defined straightforwardly making use of the relation between metrics 
as provided by a standard 4 + 1 Cauchy decomposition: 

5oo = -a 2 + 7«/3 l /3 j , goi = mP 3 , gij = 7y I C 1 - 10 ) 

where we adopt the unperturbed values for the lapse a and shift /3\ Data 
for g^t is obtained by taking a time derivative of the expressions above, 
making use of the relation between and 7^ 

-2aK lj = (d t -£p) 7ij , (1.11) 

and the harmonic coordinate condition, which in 4 + 1 form is 

(d t -(3 i d l )a = -a 2 K, (1.12) 
(d t - (3%)? = o^ l {prr n (4) r /mn - a,/) , (1.13) 

where ^Ti mn are the Christoffel symbols of the 4- metric. 

1.3.2 Evolution 

The Harmonic Decomposition with Constraint Damping 

In the following sub-section we briefly review harmonic evolution, in par- 
ticular with regards to a numerical implementation (for more details see 

4 For alternatives see jWisemanl. 120031: ISorkin and Piranl. 120031: lAnderson et all l2006fl . 
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Lindblom et al.l . 120061 ; iPretoriud . 120071 ). Harmonic coordinates are a set of 



gauge conditions that require each spacetime coordinate x a to independently 
satisfy the covariant scalar wave equation: 

V v V v x» = ~^=d v (yf^g^d^) = 0, (1.14) 



-g 

where g is the determinant of the spacetime metric g^ v H Harmonic coordi- 
nates have a long history in relativity, and are well adapted to describing 
black strings, as the metric in harmonic coordinates is regular on the hori- 
zon. More importantly, for numerical evolution, substitution of (|1.14p (and 
its first covariant gradient) into the vacuum Einstein field equations R^ u = 
yields a system of explicitly symmetric hyperbolic evolution equations for 
the metric: 

\<fi x 9^x + 3 ?X (,^)x,€ + r lx r ^ = °> ( L15 ) 

where Tf, x are the Christoffel symbols. By themselves, these equations admit 
a larger class of solutions than desired, and it is only the subset of solutions 
that satisfy what can now be considered the constraints 

C^ = g^V^x u = 0, (1.16) 

that are of physical interest. The time-derivative of these constraints can be 
related to the traditional Hamiltonian and momentum constraints discussed 
in Sec. 11.3.11 At the analytical level, Bianchi identities imply that initial 
data (the metric and its first time derivative) satisfying = as well as 
the traditional constraints, will evolve via (|1.15|) to a solution where = 
for all time, provided the boundary conditions are consistent with C M = 0. 

Numerically the situation is more complicated, as truncation error will 
generically source non-zero during evolution, and this typically grows 
exponentially (even in a convergent implementation, renderin g it difficult to 



achie ve long-time well behaved evolutions). The cure (following 



Gundlach et al. 



.2005), is to add constraint damping terms to the harmonic form of the Ein- 



stein equations (|1.15|) : 

+ k (2n^C v) - (1 + p)g^n*C x ) = 0, (1.17) 
where is a unit time-like vector, here chosen to be the vector normal to 



5 Generalized harmo nic coordinates a dd an arbitrary set of source functions H a to the right 
hand side of lITTTH lIFriedrichl. 11983) ; these can be chosen to implement different gauges, 
though here we choose = 0. 
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harmonic time t = constant surfaces, and (k, p) are the constraint damping 
parameters. Notice that the extra terms are homogeneous in C M , hence when 
Cfj_ = 0, Eqn. (|1.17|) trivially reduces to the Einstei n equations. If C u ^ , 



perturbation analysis about Minkowski spacetime (jGundlach et all 120051 ) 
reveals that for (k > 0,p > — 1), all Fourier modes of C M except a zero 
wavelength one are exponentially damped. Analytical results are not known 
for the efficacy of constraint damping in the strong-field, non-linear regime, 
though empirically it has been shown to work for generic compact binary 
systems in 4D (involving black holes and neutron stars), and the unstable 
black string in 5D. Typically, a value of k that works well is k ~ l/l, where 
i is some characteristic scale in the problem; here k is set to 1/M, where 
M is the initial mass per unit length of the unperturbed string. In all 4D 
simulations to date, p has been set to zero; for the black string, a value of 
p £ (— 1..0) was essential to damp a zero-wavelength mode growing along 
the string-like dimension w. The exact value however was not too important, 
and we chose p = —0.5. 

Regarding boundary conditions in the numerical simulation, the domain is 
periodically identified in the w direction, and at the outer boundary Dirichlet 
conditions are imposed for simplicity, with the metric fixed to the values of 
the initial data there. This latter condition is only consistent with = 
to leading order in l/rj,, where is the radius of the outer boundary; hence 
to avoid any potential "problems" that might arise from this, was chosen 
to be sufficiently large that the outer boundary is out of causal contact with 
the string horizon during the entire length of the simulation. 

Symmetries and the Cartoon Method 

Since the typical computational resources required to numerically solve a 
hyperbolic problem scale like N D , where N represents the number of mesh 
points needed to resolve a feature of interest along one dimension, and D 
is the total number of space-time dimensions, it would be impossible to 
solve for a general perturbation of a 5D (or higher) black string on con- 
te mporary computer clusters. Howeve r, as mentioned before, the analysis 
in ( Gregory and Laflamme . 1993 . 1994! ) shows that only modes along the ex- 



tra string-like direction w are unstable, whereas perturbations within the S 2 
cross sections decay exponentially. Hence it is reasonable to expect that one 
can obtain a correct picture of the final end-state by restricting to spherical 
symmetry within each w-constant slice. In other words, we consider metrics 
of the form 



ds 2 = Wg^dx^dx" + R 2 (d6 2 + sin 2 fl# 2 ) , (1.18) 
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where x^ = (x ,x ,x 2 ) = (t,r,w) (harmonic time, radial coordinate, string- 
like direction), and the three metric ^guv and areal radius R only depend 
on the 3 coordinates x^ . 

It would seem that the natural way to proceed (as similarly done with 
the initial data) is to directly discretize the field equations with a metric 
ansatz of the form (|1.18p . reducing the problem to a 2 + 1 dimensional 
simulation. However we were not able obtain long-term stable evolution 
doing so within the harmonic formulation (i.e. requiring that (t, r, w) be 
harmonic coordinates with appropriate gauge sources). Possible reasons for 
the difficulties experienced are related to spherical coordinates not being 
harmonic, and the limited set of gauge source functions we tried to ac- 
count for this fact did not yield well behaved evolutions. Potentially related 
difficulties were reported in studies of harmonic evolution in 4D spherical 
symmetry ( Sorkin and Choptuik . 201Cll ). 



A way around this problem is to consider the full 5D metric g^ u in har- 
monic Cartesian coordinates (t, x, y, z, w) as this strategy has worked well 
in 4D. The S 2 symmetry of the spacetime can then be imposed via the 
corresponding Killing vectors: 



, dy J V dx , 



dz) \dy 
dx J \dz 



e 3 M = M-l • (i-2i) 



Thus, a single (y = 0, z = 0) slice (for example) of the spacetime is sufficient 
to reconstruct the entire spacetime by the action of the Killing vectors on 
the metric: 

C*ix9w = £&9/j,v = = °- ( L22 ) 

The practical way to implemented this in the code is to discretize the given 
3D slice of the 5D metric, and then use the Killing conditions (|1.22p to 
replace derivatives orthogonal to the slice (i.e. in the y and z directions) as 
required in the harmonic evolution equations (|1.17p with derivatives tangent 
to the slice (i.e. in the x direction). For example, expanding C^g^u = one 
can solve for the ^-gradients of the metric elements as 

1 



9fiv,: 



zg»v,x - 2$ z ( P gv)x + 18 x ^g v)z . (1.23) 



This is a variant ( Pretorius . 2005bl ) of the so-called cartoon method ( Alcubierre et al 



200ll ). originally applied to axisymmetric evolution in 4D spacetime. 
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1.3.3 Apparent Horizons 

Crucial to understanding the dynamics of an unstable black string is to 
understand the behavior of its horizon. In a numerical simulation of finite 
length, at best one could recover an approximation to the actual event hori- 
zon (EH), for example by looking at the boundary of the causal past of some 
region of the spacetime at the last time step of the simulation. Furthermore, 
in spacetimes with naked singularities, as implied by the numerical solution 
for this system, an event horizon does not exist. A better local property of the 
spacetime to study is the apparent horizon (AH), defined as the outermost 
marginally outer-trapped surface (in other words the outermost surface from 
which the outward null expansion is exactly zero, and the inward null expan- 
sion is less than or equal to zero at each point on the surfa ce). Even though 



appa rent horizons are slicing dependent, the earlier study (jChoptuik et al 



2003) showed that the AH was essentially indistinguishable from the ap- 



proximate EH everywhere, modulo at late times because of the ambiguity 
in defining the region of spacetime whose causal past defines the exterior of 
the EH. We thus focus on the AH, which incidentally is also essential for 
the excision technique used to remove the geometric singularity inside the 
black string from the computational domain. 

For more informatio n on how AHs ar e defined and searched for in a nu- 
merical evolution, see (|Thornburg| . l2007h : we use a flow method in this code. 



1.3-4 Evolution Code 

The basics of the numerical evolution code employed can be briefly summa- 
rized as follows. The harmonic equations with constraint damping ()1.17p . 
reduced to first order in time by introducing the auxiliary variables g^ u = 
dt.Q/uj, are discretiz e d usin g 4 th order finite difference methods with AMR 



( Berger and Oliger . 1984 ). as implemented in the PAMR/AMRD libraries 



(which also handle parallelization via MPI). Time integration is via 4 th or- 
der Runge-Kutta. Standard, centered spatial difference operators are used 
in the interior of the grid, while at the inner excision boundary (chosen to 
match the shape of the AH, but a fraction between 10 and 50% smaller in 
radius), centered difference operators are replaced with sideways operators 
as appropriate. A 6 th order Kreiss-Olig er style dissipation filt er is used to 
control high-frequency truncation error (jCalabrese et all l2004h . 



http : //laplace .physics .ubc . ca/Group/Sof tware . html 
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1.4 Evolution of an Unstable, 5D Black String 

With the previously described implementation, we now concentrate on de- 
scribing the results from the evolution of an unstable black string. We con- 
centrate on a single case where the unperturbed string has periodicity length 
20M. Such an identification satisfies L > L c , where L c is the critical wave- 
length beyond which modes become unstable, yet only allows a single un- 
stable mode initially. To study this case, we adopt a computational domain 
(r,w) G ([0, 320M] x [0,20M]) (in the following for convenience we relabel 
the x Cartesian coordinate as r, though of course x = r on the y = z = 
slice of the spacetime). We adopt an initial grid that uniformly covers the 
entire domain with (N r ,N w ) = (1025,9) points. As the evolution proceeds 
and finer structure arises, the AMR algorithm introduces additional higher 
resolution grids where needed, based on a specified maximum truncation 
error tolerance. 

A typical "low" resolution run has a coarsest radial mesh spacing of ~ 
0.3/M, which is the resolution at the outer boundary. Additional levels are 
added by the AMR algorithm (each refined with a 2 : 1 ratio, in both r and 
w, of the parent level) to resolve the region of spacetime near the AH of 
the string. Initially 3 additional levels are added, growing to 16 additional 
levels by the time the simulation was terminated at t ~ 229. 0M (due to the 
prohibitive computational costs that would have been required to continue). 
"Medium" and "high" resolutions, for convergence and error estimates, were 
specified by decreasing the maximum allowed (estimated) truncation error 
by a factor of 8 and 64 respectively. Due to lengthier run times with increased 
resolution, the medium and higher resolution cases were not run for as long. 
Specifically, the longest medium (high) resolution case was terminated at 
t ~ 226.8M (t ~ 217.9M); the high resolution run required around 100,000 
CPU hours on the woodhen cluster at Princeton University 0. This translated 
into roughly 2 months of essentially continuous running on 100 processors. 



1.1^.1 Apparent Horizon Dynamics 

To help understand the dynamics of the system, we monitor several relevant 
quantities, in particular, the apparent horizon radius RAH(t,w), the total 
horizon area A(t), its intrinsic geometry via a flat-space embedding diagram, 
as well as two spacetime curvature invariants I = R^ UX ^R^ UX ^ and J = 
R /11/x ^R x ^ ricr Rrja^ evaluated on the horizon. For the curvature invariants, 

7 http : / /www . pr inceton . edu/ resear chcomput ing/comput at ional -hardware/ 
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Figure 1.1 Apparent horizon area as a function of time for the perturbed 
black string. What is shown is the area normalized to the initial area, and 
from simulations with 3 different resolutions (see the text for a description 
of the resolutions) . 



we find it useful to rescale them as, 

K = IR\ H /Yl, S = 27 (l2J 2 r 3 - l) + 1 , (1.24) 



as such rescaling yields K = S = 6 for an S s Tangherlini black hole while 
K = S = 1 for a uniform black string. 

Fig. Il.ll shows the total apparent horizon area A as a function of time for 
the evolution of our perturbed black string. As expected for a reasonable 
approximation to the event horizon, the area is non-decreasing with time. 
More interestingly, we note that at the end of the simulation corresponding 
to the lowest resolution run (the one that ran the farthest in time) the total 
area is A = (1.369 ± 0.005)A)J! 

where Aq is the initial area; this value 
essentially saturates the value of 1.374j4o that an exact 5D black hole of the 
same total mass would have. 

More insight into the dynamics of the horizon can be garnered by observ- 
ing the evolution of its intrinsic geometry with time. Fig. 11.21 shows several 
snapshots of embedding diagrams of the AH, from the medium resolution 
simulation run. As can be seen in the figure, the string initially evolves to 
a configuration resembling a hyper sp herical black hole conn ected by a thin 
string segment, as reported earlier in (jChoptuik et all . 1200.4 where the sim- 
ulation ended at roughly t ~ 164M). However, as the evolution proceeds 
it is apparent that the string segments are themselves unstable, and this 
pattern repeats in a self-similar manner to ever smaller scales. 

Though the intrinsic geometry on the horizon resembles spheres connected 
by strings, this by itself does not imply the local spacetime geometry is simi- 

8 The error in the area was estimated from convergence at the latest time data was available 
from all simulations. 
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lar to either 5D spherical black holes or black strings, respectively. However, 
evaluation of the curvature invariants (jl,24p on the horizon gives further in- 
sight into this question, and this is shown in Fig. ll.3t taken from the last time 
step of the medium resolution run. As seen in the figure, near spherical-like 
sections the normalized invariants approach 6, in contrast to the string-like 
sections where they are close to 1; thus, at least based on these two invariant 
indicators, the near horizon geometry indeed resembles that of the solutions 
suggested by the shape of the AH. 

1.4-2 Interpretation of the horizon dynamics 

The shape of the AH in the embedding diagram, and that the invariants 
are tending to the limits associated with pure black strings or black holes 
at corresponding locations on the AH, suggests it is reasonable to describe 
the local geometry as being similar to a sequence of black holes connected 
by black strings. This also strongly suggests that satellite formation will 
continue in a self-similar cascade, as each string segment locally resembles a 
rather uniform black string, and is sufficiently thin and long to be unstable. 
Note that even if at some point in the cascade thick segments were to form, 
this would generically not be a stable configuration, as (except possibly 
in an exactly symmetric situation) the satellites will have some non-zero 
tt>-velocity; hence they would eventually merge, effectively lengthening the 
segments connecting the remaining satellites. 

With this interpretation, we can understand key features of the AH dy- 
namics by sequencing its time evolution into generations of satellite forma- 
tion. We assign a time t{ to the onset of a new generation when the local 
instability has reached an "observable" level, defined here (somewhat arbi- 
trarily) as the time when a nascent spherical region reaches an areal radius 
1.5 times the surrounding string radius. For each generation we measure the 
number n s of satellite black holes that form per string-segment, their radii 
Rah,i as well as the corresponding string segment's radii R S) i and lengths 
L Sj j. These features are summarized in Table [TTTl where error bars, where 
appropriate, come from convergence calculations^. Only the first generation 
is in a sense non-generic, as, by construction, only one unstable mode is 
initially allowed to develop. Generations after the first one display string 
segments where their ratio L 8 i/R a i ~ 100 is not only above the critical one 
(~ 7.2), but can in principle accommodate several unstable modes (though 

9 The global truncation error tends to grow with time; also, as mentioned, for the 3rd 

generation we only ran the low and medium resolution cases far enough to uncover it, hence 
(assuming convergence) the much larger errors there. Only the low resolution run was 
continued for long enough to reveal a 4th generation, thus the lack of error-estimates there. 
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Figure 1.2 Embedding diagram of the apparent horizon at several instances 
in the evolution of the perturbed black string, from the medium resolution 
run. R is areal radius, and the embedding coordinate Z is defined so that 
the proper length of the horizon in the space-time w direction (for a fixed 
t, 8, (f>) is exactly equal to the Euclidean length of R(Z) in the above figure. 
For visual aid copies of the diagrams reflected about R = have also been 
drawn in. In the top panel, the thin black lines are from the initial time; 
in the top two panels the thick black lines are the last time from the time- 
segment depicted in the corresponding panel, whereas in the bottom panel 
the last time-snapshot has been drawn with a thin black line. Note that the 
vertical and horizon axis scales have been changed in the bottom panel to 
better show the satellites that have formed at late times. The computational 
domain is periodic in w with period 5w = 20M; at the initial (final stage 
of the simulation) SZ = 20M, (SZ = 27.2M)— see also Fig. O 



see the discussion in Sec. 11.5. 1|) . One would expect the mode closest to the 
maximum in the dispersion relation to dominate in each segment, and thus 
qualitatively similar dynamics to unfold from one generation to the next. In 
particular, we see that the time scale for the development of the n th gen- 
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Figure 1.3 Curvature invariants evaluated on the apparent horizon at the 
last time of the medium resolution simulation (as depicted in Fig. II. 2[) . 
Also, for visual aid, the areal radius of the apparent horizon is shown. The 
invariant K evaluates to 1 for an exact black string, and 6 for an exact 
spherical black hole; similarly for S (|1.24l) . 
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R s ,i/M 


Rahj/M 


Ls.i/ Rs,i 


1 


118.1 ±0.5 


1 


2.00 


4.09 ± 0.5% 


10.0 


2 


203.1 ±0.5 


1 


0.148 ± 1% 


0.63 ± 2% 


105 ± 1% 


3 


223 ±2 


> 1 


0.05 ± 20% 


0.1-0.2 


« 10 2 


4 


« 227 


> 1 


« 0.02 


? 


« 10 2 



Table 1.1 Properties of the black string apparent horizon (see the text for a 

discussion). 



eration of the instability is essentially the same along each string segment, 
as are the radii of the corresponding strings and spheres that form. How- 
ever, beyond the second generation we find variation in the number n s of 
satellites that form as a function of resolution (hence the > 1 in the table). 
This suggests that exactly which mode dominates depends sensitively on 
the initial conditions, sufficiently so that a small perturbation, here coming 
from numerical truncation erroj 10 !. can change the quantitative details of the 
late time dynamics. This is not unexpected, as from the Gregory-Laflamme 
dispersion relation, except for the mode exactly at the maximum, there are 
two unstable modes with the same growth rate. 



In slight abuse of terminology, as of course truncation error does not in general correspond to 
any physical perturbation. 
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Extrapolation to the end-state 

The above observation of solution properties allow us to extrapolate to the 
end-state of the instability, estimating the time when this will occur, and 
the structure of the horizon just prior to it. We first calculate when the 
self-similar cascade will end, namely, the time when the connecting string 
segments reach zero radius, and then estimate the fractal dimension of the 
AH geometry just prior to this end. 

The time when the first generation satellite appears is controlled by the 
perturbation imparted by the initial data, which here is Tq/M ~ 118. Subse- 
quent generations, however should represent the generic development of the 
instability. From the data in the table, the time of growth of the first insta- 
bility beyond the one sourced by the initial data is T\/M ~ 80. Beyond that, 
with the caveats that we have a small number of data points and poor con- 
trol over the error at late times, the data suggests each subsequent instability 
unfolds on a time-scale X ~ 1/4 times that of the preceding one. Again this 
is to be expected if the instability is qualitatively like the Gregory-Laflamme 
instability of the exact black string, where the time scale is proportional to 
the string radius. Then, the total time AT for the end-state to be reached 
is roughly 

OO rp 

AT^T + J2TiX l = T + — 1 —. (1.25) 

For this case then, AT/M ~ 231. At this tim^] all string segments will 
reach zero radius. Since the Kretchman curvature invariant (|1.24p just out- 
side a string segment of radius r is proportional to r -4 , and recalling that 
here (harmonic) time is regular everywhere from some distance inside the 
AH outwards and corresponds to the time measured by stationary asymp- 
totic observers, this indicates formation of a naked singularity, and thus a 
violation of cosmic censorship. Furthermore, it is generic in the sense that no 
fine-tuning of the initial data is required (i.e., any length- wise perturbation 
exceeding the critical wavelength will do). 

To give further evidence that the data supports the above conclusion, in 
Fig. 11.41 we show the time evolution of the radius of several representative 
cross-sections of the AH, in logarithmic coordinates where time has been 
shifted by the estimated time of naked singularity formation. That the ra- 
dius of string-like segments decrease linearly in shifted-logarithmic time is 
consistent with a self-similar scaling to zero radius at the corresponding fi- 

11 The exact value of which could be expected to vary slightly along the length of the horizon 
with our chosen time coordinate and initial perturbation, though one should be able to define 
a time-slicing where this time is exactly the same everywhere along the string. 
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nite asymptotic time. Another intriguing aspect of this self-similar scaling is 
that it is qualitatively the same as that observed in the approach to pinch-off 
of Raleigh-Pla teau unstable fluid streams. In the fluid scaling solu- 

tion is known (jEggersl . Il993l : iMivamotol . l201(lh , and the radius R of the fluid 



column decreases linearly with time t to fluid breakup t = t c : 

R<x(t c -t). (1.26) 

Fig. 11.41 shows that this, to good approximation, also describes shrinking 
regions of the black string. 

A further consequence of the self-similar nature of the instability is that 
the AH shape will develop a fractal structure prior to pinch-off. This implies 
that for a fixed angular cross section (8, (ft = const.) of the horizon, the proper 
length L p (t) (within the periodically identified domain) will grow with each 
subsequent generation, diverging at a rate related to the Hausdorff dimension 
d of the end-state shape. Assuming a scaling relation of the form (jl.26p . 
that the additional length in each new generation is caused by every string 
segment developing (on average) the same number n s of satellites, each 
with radii Rahj following the scaling suggested in Table II. 1^ and that n s x 
Rah J /L s ,i <C 1, one can show that the following growth of L p (t) is expected 

L p (t)oc(t c -t)^- d l (1.27) 

Figure [T31 shows a plot of L p (t) on a logarithmic scale. From the measured 
slope and above relationship, we find d ~ 1.05, which as expected is greater 
than a value of 1 (that would correspond to a non-fractal curve) but only 
slightly, as the fractal structure is obtained by repeatedly replacing a rela- 
tively long straight line with a similar length line plus a small semi-circular 
protrusion. 



1.5 Speculations and Open Questions 

In this chapter we described the evolution of an unstable black string in 
5D, at least as far as was numerically feasible given finite computational 
resources, and extrapolated the observed behavior to describe the nature 
of the end-state of the instability. This has provided a first insight into 
the fascinating dynamics of unstable horizons in higher dimensions, yet has 
raised many questions and opened up avenues for future work. In this section 
we conclude by discussing in Sec. 11,5. Tl whether the Gregory-Laflamme linear 
perturbative analysis can be used to quantitatively understand the dynamics 
of subsequent generations of the string instability, in Sec. 1 1 . 5 . 21 mention when 



1.5 Speculations and Open Questions 



19 




Figure 1.4 Logarithm of the areal radius vs. logarithm of time for select 
points on the apparent horizon, from the same (medium) resolution run 
depicted in Fig. 11.21 We have shifted the time axis assuming self-similar 
behavior; the putative naked singularity forms at asymptotic time t/M rs 
231. The coordinates at w = 15,5 and 4.06 correspond to the maxima 
of the areal radii of the first and second generation satellites, and one of 
the third generation satellites at the time the simulation stopped (the w 
positions of the satellites do slightly evolve with time). The value w — 6.5 
is a representative slice that is in the middle of a piece of the horizon that 
remains string-like throughout the evolution. For visual aid in comparison 
with the analogue fluid pinch-off scaling solution (|1.26p . a line with slope 
equal to -1 has been added (thin brown line). 




-5 -4 -3 -2 -ln(231-t/M) 



Figure 1.5 The proper length L p (t) vs. time of the apparent horizon in 
logarithmic coordinates, along a fixed azimuthal cross-section (9 = <j> = 
constant in spherical coordinates) and over one period of the spacetime in 
the w direction. The time coordinate is scaled as in Fig ll.4l The dotted line 
is a visual linear fit to the late-time dynamics, and has a slope of ~ 0.048. 



20 



Final State of Gregory-Laflamme Instability 



quantum corrections are expected to become important as the string thins, 
and end with a list of open questions for future directions in Sec. 11.5,31 

1.5.1 Mode behavior 

To gain more insight into what determines the observed dynamical behavior, 
we analyze the different modes in a string segment at the late stages of the 
first generation shown in Table 11.11 corresponding to the early stages of 
the growth of the second generation. In particular, we look at RAH(t,w) in 
t 6 (175M, 205M) when only one spherical region is clearly distinguishable, 
and, in our coordinates this region has extent inwG [10, 20M]. The string 
section corresponds to w E [0, 10M], and we decompose R^H(t,w) in this 
latter domain via the following expansion 

R AH (t, w)=c - Eg lCI sin(/™) , (1.28) 

and extract the coefficients q up to I = 6. Figure 11.61 illustrates the val- 
ues q/co within this time frame, and we only plot the odd values of I as 
{c2,c^,cq} are about two orders of magnitude smaller. This is consistent 
with the observed symmetry in the development of the single second gen- 
eration hypersphere within this segment, as it forms exactly in the middle 
of the string. Interestingly, of the odd modes, only I = 1 displays growing 
behavior even though the length/radius ratio of the string at this time would 
allow for several modes to be unstable according to the Gregory-Laflamme 
instability criteria. This may not be too surprising however, as during this 
stage the solution is highly dynamical, shrinking in radius at a rate compa- 
rable to that at which putative Gregory-Laflamme modes could grow. This 
suggests that the linear analysis perturbing about the static black string 
"background" is not applicable here, and can only capture the qualitative 
nature of subsequent generations of instability. Note though that this is not 
in tension with the observation provided by the invariants K and S, as the 
string can shrink in a time-dependent fashion and still maintain K ~ S ~ 1. 

1.5.2 Dynamics Beyond the Classical Regime 

As discussed, our classical description of the system reveals cascading self- 
similar behavior of the horizon, with regions evolving to ever shrinking cross- 
sectional radius. Such regions asymptote (in finite time) to a zero-mass naked 
singularity. Obviously for small enough regions, a classical description is no 
longer applicable and quantum phenomena must be taken into account. This 
is typically expected to happen when the cross-sectional radius r approaches 
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Figure 1.6 Vales of expansion coefficients for the string segment corre- 
sponding to the first generation. Even modes (except Co) are a few orders 
of magnitude smaller than odd modes. 



some small length I scale (e.g. the Planck length l p or the string length l s ), 
though exactly when depends on the details of the leading order corrections 
to the field equations. For example, stringy-corrections predict this to be of 
the form 

Ruu = fRua^Ru^ + • • • • (1-29) 

Since each term in (jl.29p has dimension 1 /length 2 , a naive estimate of when 
the RHS becomes important is when its magnitude becomes of order 1/r 2 , 
as r is the relevant length scale in the vicinity of a shrinking region of 
the horizon. Since R^a^ ~ V 7 " 2 there, this suggests corrections to general 
relativity are important when £ 2 /r 4 ~ 1/r 2 ; i.e. if t = l p , when the radius 
shrinks to the Planck length. 

Another effect that could in principle alter the classical picture before the 
Planck scale is reached is if the timescale of Hawking evaporation tjj of the 
string segment becomes smaller than the Gregory-Laflamme timescale tql- 
However, as the following shows, th only becomes smaller than tql beyond 
the Planck scale. 

Assuming that a 4D description of evaporation at a w = const, cross 
section of a string-like region gives a decent approximation to an evaporating 
string, we can estimate the timescale of Hawking evaporation and compare 
it to the timescale for development of the next generation within the self- 
similar cascade. From the simulation results, we have roughly that 



t gl ~ WOMiG/c 3 , 



(1.30) 
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while that for Hawking evaporation of a 4D black hole of mass M, is 

Thus, for these two times to become comparable 

Mi~ 10-M— = 10- x Mp, (1.32) 

where Mp is the Planck mass. This suggests that if the length scale t in (|1.29p 
is of order the Planck length or larger, then higher curvature corrections will 
generically be more important in describing the leading order alterations to 
the classical description of the instability than Hawking evaporation. 



1.5.3 Future work 

We conclude this chapter by listing some open questions and directions for 
future work, though this is by no means an exhaustive list. 

First, consider a broader class of initial perturbations of the 5D black 
string where other unstable modes, and possibly linear combinations of 
modes with similar growth rates, are exited. We do not expect anything 
to change in a radical fashion, as the qualitative nature of the linear growth 
of all these modes is similar. Nevertheless, a quantitative analysis of such 
scenarios would tell if this expectation bares out, and allow one to under- 
stand how (if at all) the initial perturbation affects the unfolding structure 
of subsequent generations of the instability. On a similar vein, exploring how 
adding angular momentum to the string changes the picture would be in- 



m 



in 


Dias et al. 


2010, 


Marolf and Palmer 



20041 it was suggested that rotation could induce super-radiant and gyrating 
instabilities; this should make for much richer dynamics in the approach to 
the end-state. Allowing for angular momentum would require less symmetry, 
making for more expensive numerical evolution. 

Second, extract more details about the spacetime, such as the evolution 



of the generators of the horizon (as in iGarfinkle et~ail 120051 ). and the gravi- 



tational waves that are emitted as the instability unfolds. Perhaps the most 
profound question raised by these results is the striking qualitative similarity 
between the non-linear evolution of the instability and the Rayleigh-Plateau 
instability. Is this a coincidence, or the consequence of a deeper relationship 
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between Einstein and Navier-Stokes than already suggested by the mem- 
brane paradigm, black folds and other perturbative descriptions of horizon 
dynamics? What would be useful in trying to understand this is to identify 
and extract geometric characteristics of the dynamical horizon that would 
map to effective fluid properties. Such a map could also prove useful in 
taking what is known about the Rayleigh-Plateau instability to learn more 
abou t Gregory-La flamme; for example, translating the Eggers scaling solu- 
tion ( Eggers . 19931 ) to an approximate solution for a thinning string segment. 

Third, consider non- vacuum spacetimes; in particular charged black strings 
would b e interesting as charges modify the onset of instabilities in t he system 
(for e.g. Gregory and Lanamme . 1994 ; Sarbach and Lehner , 20051 ). 

Fourth, study where the spacetime is asymptotically flat in all 

spatial directions (the current study imposed periodicity in the fifth dimen- 
sion). This would be required to provide an example of cosmic censorship 
violation in the 5D asymptotically flat case. However, it is natural to expect 
the same behavior as seen here based on the following putative scenario. 
Begin with a highly distorted S 3 horizon, namely one which is long and 
thin ("cigar shaped") so that near the center of the horizon it locally re- 
sembles an S 2 x R black string. One may expect that this horizon would 
ring-down via gravitational wave emission to a uniform S* 3 horizon; however, 
if the dynamical time-scale for the ring-down (which will be proportional to 
the light-crossing time in the prolate direction, so could be made arbitrarily 
long) is much longer that the Gregory-Laflamme instability of the central 
region, then the latter should take over first, resulting in a pinch-off. 

Fifth, explore black strings in higher dimen sional spacetimes, particularly 
to investigate the conjecture of (jSorkinl . l2004j ) that a qualitatively different 
end-state is expected for D > 13. If one continues to impose SO(D — 2) 
symmetry the problem can still be expressed to depend only on (t, r, w), 
making it tractable with current computational resources. However, as the 
fields decay as {r D ~~ 3 ,r D ~ 2 } moving away from black hole, black string so- 
lutions respectively, the resolution and/or order of the numerical scheme 
would need to be higher than the one used in this work to obtain results of 
comparable accuracy. 

Sixth, explore additional black objects subject to Gregory-Laflamme-like 
instabilities. For instance, rapid ly rotating black holes have been shown to be 
subject to a similar instability (jEmparan and Mversl . I2003T). Recent numeri- 
cal w ork presented the first exploration of such systems (jShibata and Yoshind . 
20ld ). though there radiation of angular momentum stabilizes the black 



holes considered. However, the growth rate of the instability in the cases 
studied was rather mild, and it is quite likely that choosing a more extreme 
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scenario (where the timescale of the unstable modes would be shorter than 
the dynamical time of gravitational radiation) would give rise to a pinch-off 
similar to the one studied here. 

Finally, evolve unstable black strings (and other black objects) within 
asymptotically Anti de Sitter (AdS) spacetime, and explore the consequences 
within the context of the AdS/Conformal Field Theory (CFT) dualities of 
string theory. This is partly related to the point mentioned above about fur- 
ther investigating the intriguing connections between gravity and fluids, as 
certain states within a CFT will admit a hydrodynamic description. Though 
regardless, understanding the CFT-duals to unstable black hole spacetimes 
is interesting in its own right, and could have implications to the more re- 
cent applications of AdS/CFT to model certain condensed matter and high 
energy particle physics systems. 
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